function [decompOutput, decompVars, decompShocks] = ...
    perform_spectral_decomp(StructureM_, Structureoo_)

% Frequency Decomposition
% Code loads information from DYNARE output and uses code by Alejandro
% Justiniano to perform frequency decomposition

% Set path to dynare as some dynare functions are called. Works with
% version 4.2.0
% See also detailed description in vardecom_spec.m

% CALLS:
% - find_dynare_transition_eq_matrices
%       - kalman_transition_matrix
% - build_selection_matrix
% - time_to_frequency
% - vardecom_spec
%       - lyapunov_symm

% NOTE: - MAKE VARIOUS SPECIFICATIONS IN ALL SECTIONS

% Christoph Goertz
% 17/04/2012



%% Inputs needed (from DYNARE): 
% GG,HH,SDX,ZMAT correspond to the solution of the DSGE Model
% TRANSITION Equation is
% $$ S_{t}= G1 * S_{t-1} + HH * e_{t} $$ where 
% $$ V( e_{t} ) = (SDX')*SDX $$ 
% Observation Equation is 
% $$ Y_{t} = zmat*( S{t} + C ) $$ and C is not used here. 
%
% FREQ:      Vector of frequencies
% EDGEW:     Distance of edge to border (2 x 1 vector)
% STLOC:     Location of rows of the state vector wish to analyze
% LEVINDIC:  Indicator, ==1 will compute the level by applying difference
% filter


%%                                Output
% Dimension NF=# frequencies, NST=# states; NLEV=# states transformed to
%           levels
%
%1. VCVZ & VCVST    Usual variances for observables and selected states
%   VCVZ cannot be generated as necessary function is not available. See
%   vardecom_spec.m (section 2)
%
%2. SPEC: Spectrum matrcies at each frequency
% -------------------------------------------
% SPECZ           [NZ NZ NF]    For observables
% SPECST          [NST NST NF]  For untransformed states
% SPECLEV         [NLEV NF]     If states transformed to levels,
%                 Diagonal elements only
%
% 3. VFB: total variance for this frequency band
% ------------------------------------------------
% VFBZ          [NZ 1] (observables)
% VFBST         [NST 1](unstranformed states)
% VBFLEV        [NLEV 1] Total variance for level states
%
% 4. SPECDEC    Variance decompsition for that frequency
%               Diagonal elements ONLY
% -------------------------------------------------------
% SPECDECZ      [NZ NX NF]
% SPECDECST     [NST NX NF]
% SPECDECLEV    [NLEV NX NF]
%
% 5.VFBDECL total Variance decompositon for this frequency band
% Ratio of the total variance in SPECDEC to the total variance in
% VFB
% ------------------------------------------------------------
% VFBDECZ      [NZ NX]
% VFBDECST     [NST NX]
% VFBDECLEV    [NLEV NX]


%% LOAD AND SPECIFY INFORMATION TO GENERATE INPUTS
% Assign information from estimation output
M_ = StructureM_;
oo_ = Structureoo_;

% Specify variables of interest and observables variables to find
% corresponding index number in dr.
 observablesStr = [
           'dy        ';
           'dc        ';
           'dinve     ';
           'labobs    ';
           'dw        ';
           'robs      ';
           'pinfobsc  ';
           'pinfobsi  ';
           'dtfp      ';
           ];
 variablesStr = [
           'dy        ';
           'dc        ';
           'dinve     ';
           %'lab       ';
           'dw        ';
           %'r         ';
           %'pinfc     ';
           %'pinfi     ';
           %'spreadc   ';
           %'spreadi   ';
           'dtfp      ';
           ];
% varStr = [observablesStr; variablesStr];
% % Find corresponding index in M_.endo_names
% for iVarStr = 1:size(varStr,1)
%     for iEndo_names = 1:size(M_.endo_names,1)
%        strInd = strcmp(strcat(varStr(iVarStr,:)),...
%            strcat(M_.endo_names(iEndo_names,:)));
%         if strInd == 1
%             varIx(iVarStr) = iEndo_names;
%         end
%     end  
% end
% if size(varStr,1) ~= size(varIx)
%     warning('Problem in variable specification section');
% end
% % Since indices are different in the dr structure, find corresponding
% % indices
% varIx_dr = oo_.dr.order_var(varIx);
% 
% % Specify indices of the observables (indices are associated with order in
% % observablesStr)
% ixObservables = varIx_dr(1:size(observablesStr,1));
% % Specify indices of (non-observables) variables to be decomposed (indices
% % are associated with order in variablesStr)
% ixVariables = varIx_dr(size(observablesStr,1)+1:size(observablesStr,1)+size(variablesStr,1) );


% INDICES NEED TO BE SPECIFIED ACCORDING TO ORDER IN oo_.SmoothedVariables!

% dy dc dinve labobs dw robs pinfobsc pinfobsi dtfp
%ixObservables = [14, 15, 16, 19, 18, 20, 21, 22, 51 ];
ixObservables = [8, 9, 10, 13, 12, 14, 15, 16, 20 ];
% dy dc dinve dw dtfp
%ixVariables = [14, 15, 16, 18, 51];
ixVariables = [8, 9, 10, 12, 20];

% Indicate which variables in variablesStr will be decomposed in levels
% (not an option for observables)
inLevels = [1 1 1 1 1];
%inLevels = zeros(10,1);
clear ix strInd iEndo_names iVarStr

%% GENERATE INPUTS FOR SPECTRAL DECOMP FUNCTION

% 1. Find matrices of the transition matrix
[A,B] = find_dynare_transition_eq_matrices(oo_, M_);
G1 = A; % matrix of predetermined variables effects in linear solution
HH = B; % matrix of shocks effects in linear solution


% 2. SDX is a matrix of zeros with of the shocks' standard deviations on
% the diagonal, as it needs to fulfill  V( e_{t} ) = (SDX')*SDX (where
% V(e_{t}) is the variance of the exogenous variables).
SDX = sqrt(M_.Sigma_e);


% 3. Use dynare function to build selection matrix of observation equation
mf = ixObservables;   % p*1 vector, indices for the observed variables
m = length(G1);       % number of endogenous variables
p = length(mf);
Z = build_selection_matrix(mf,m,p);
zmat = Z;


% 4. Time needs to be transformed in frequencies
freq = time_to_frequency(6,32,500);         


% 5. Specify distance of edge to the lowest and highest frequency 
% respectively
edgew = [freq(1)-freq(2); freq(end-1)-freq(end)]/2;     


% 6. Location of rows of the state vector wish to analyze
stloc = ixVariables;


% 7. Indicate which variables will be decomposed in levels
levindic = inLevels; % needs to be same length as stloc.

clear A B ys ExoValue exoValue iexoNames mf m p Z

%% PERFORM DECOMPOSITION AND SPECIFY OUTPUT
% Only the first six inputs are needed when decomp is performed only for
% observables.
out = vardecom_spec(G1,HH,SDX,zmat,freq,edgew,stloc,levindic);

% Specify decomposition output according to variables chosen above
displayColumns = [1:7,10:12,14:15,28:29,31:32]; % Display only for nonzero shocks
decompObservables = out.vfbdecz(:,displayColumns); % Decomp of observables 
                                                   % (not in levels)
decompLevels = out.vfbdeclev(:,displayColumns); % Decomp of variables 
                                                % in levels
% Replace some observables by decomp in levels:
decompOutput = [decompLevels(1:3,:);decompObservables(4,:);...
    decompLevels(4,:);decompObservables(6:8,:);decompLevels(5,:)];

% Corresponding variable names:
decompVars = [variablesStr(1:3,:);observablesStr(4,:);...
    variablesStr(4,:);observablesStr(6:8,:);variablesStr(5,:)];
% Corresponding shock names: 
decompShocks = M_.exo_names(displayColumns,:);

% Corresponding variable names: observablesStr, variablesStr
% Corresponding shock names: M_.exo_names(displayColumns,:)



